Parametric resonance of capillary waves at the interface between two immiscible Bose-Einstein 

condensates 
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We study parametric resonance of capillary waves on the interface between two immiscible Bose-Einstein 
condensates pushed towards each other by an oscillating force. Guided by analytical models, we solve numeri- 
cally the coupled Gross-Pitaevskii equations for two-component Bose-Einstein condensate at zero temperature. 
We show that, at moderate amplitudes of the driving force, the instability is stabilized due to non-linear modifi- 
cations of the oscillation frequency. When the amplitude of the driving force is large enough, we observe detach- 
ment of droplets from the Bose-Einstein condensates, resulting in generation of quantum vortices (skyrmions). 
We analytically investigate the vortex dynamics, and conditions of quantized vortex generation. 



I. INTRODUCTION 

Numerous recent studies of hydrodynamic phenomena in 
quantum media have demonstrated remarkable interplay be- 
tween quasi-classical hydrodynamics and purely quantum ef- 
fects such as quantum solitons and vortices, e.g. see |Q]-|6l]. As 
an example of such interplay, we can mention works on shock 
fronts in Bose-Einstein condensates (BECs), non-linear optics 
and quantum plasmas |7|-|9|]: the quantum shocks propagate as 
soliton trains due to Bohm-de Broglie dispersion instead of the 
monotonic transition between low-density and high-density 
gases (fluids, plasmas) in a classical shock. Development of 
quantum solitons has been also obtained experimentally and 
theoretically in the process of dynamical quantum interpene- 
tration of both miscible and immiscible BECs I2I ITol [TTfl . 

Another important example of the interplay concerns gen- 
eration of quantized vortices. Traditionally, quantum vortices 
in BECs are produced by rotating the trap [12], by stirring 
BECs by moving potentials B13I1 . by coherent transfer of the 
orbital angular momentum of photons to atoms 1U4II . by adia- 
batic phase imprinting 01511 . or by modulational instability of 
solitons lfl6[|r7ll . Recently, there has been also much interest 
in multi-component BECs with a well-distinguished interface, 
which allow the possibility of vortex generation by means 
of quasi-classical hydrodynamic instabilities: the Kelvin- 
Helmholtz instability, the Rayleigh-Taylor (RT) instability, 
the capillary instability, etc. 1141461 1 18142 ill . With the help of 
the instabilities, one can produce rather complicated quantum 
vortex structures in BECs like skyrmions, for which the intrin- 
sically empty vortex core in one BEC-component is filled by 
the other component lfTlll2lll . The system of a two-component 
BEC can be realized experimentally by blocking the spin- 
recombination process 1 1 , — 1) + 1 1 , 1 ) — > 1 1 , 0) + 1 1 , 0) by the 
quadratic Zeeman effect, as has been discussed in iflal . In 
the case of magnetic gradient pushing the BEC components 
to each other, they tend to reduce the potential energy of the 
system by exchanging places, which typically happens in the 
form of a multidimensional flow due to the RT instability. In 
that sense, the magnetic field gradient plays the role similar to 
a gravitational field for the classical RT instability. As an alter- 
native to the RT instability, Ref. Hill has found the possibility 
of reducing the excess of potential energy by ID dynamical 
quantum interpenetration; still, it has been also shown that the 



RT instability dominates for interface perturbations of a suffi- 
ciently large wavelength. 



Linear stability analysis based on the variational principle 
1I20I1 has demonstrated that the same experimental configura- 
tion may also reproduce other quasi-hydrodynamic phenom- 
ena in BECs such as the Richtmyer-Meshkov instability and 
the parametric instability by using time-dependent magnetic 
field, which corresponds to time-dependent effective "grav- 
ity". We focus on the parametric resonance, which is also 
known as the parametric or Faraday instability. We point out 
that interface dynamics driven by time-dependent "gravity" 
is intensively studied within the classical fluid (gas) mechan- 
ics, e.g. see recent numerical simulations of the Richtmyer- 
Meshkov instability in gases I22I - I24I1 . laboratory experiments 
on fluids with oscillating effective "gravity" 11251 uql. w orks 
on inertial confined fusion l27l l28ll and combustion ll29l - [32ll . 
In the classical hydrodynamics, the parametric instability is 
typically excited at the interface between light and heavy 
gases (fluids) with the "gravity" produced by direct vibrations 
of the experimental set-up or by acoustic waves (the later is 
especially typical for combustion systems). In particular, ex- 
perimental and numerical studies of the parametric instability 
in combustion encountered powerful turbulence produced in 
the flow because of the instability JHHHH]. 



In the present paper we study development of the paramet- 
ric instability at the interface of two immiscible BEC com- 
ponents pushed towards each other by an oscillating force. 
The instability arises due to the parametric resonance pump- 
ing quantum capillary waves at the interface. We show that 
in the present configuration the instability does not lead to 
turbulence. At moderate amplitudes of the driving force the 
instability is stabilized at the nonlinear stage due to modi- 
fications of the doubled oscillation frequency in comparison 
with frequency of the driving force. When the amplitude of 
the driving force is large enough, we observe detachment of 
droplets from BEC components and generation of quantum 
vortices (skyrmions). We discuss properties and dynamics of 
the skyrmions. 



H. ANALYTICAL MODEL FOR THE PARAMETRIC 
RESONANCE 

Two immiscible BECs are separated in space in the ground 
state of the system. To investigate physics of the dynami- 
cal system it is sufficient to study a symmetric case, where 
the components' atoms have equal masses, and equal intra- 
component interaction parameters g = g\\ = g22, gn > g, 
where gjj = 47th 2 aij/m and a, ; - are the scattering lengths for 
collisions between atoms of the z'-th and j-th components. 
The scattering lengths are calculated depending on the geom- 
etry of system confinement J35ll : the dimensionless parame- 
ter 7 = gn/g — 1 > specifies strength of mutual repulsion 
of the components. We study a ribbon-shaped 2D geometry 
where the BEC components are initially placed in the domains 
y < and y > 0, respectively. The BEC is tightly confined 
in the z direction; in the y direction it is trapped so that the 
Thomas-Fermi (TF) approximation for each component holds 
well 111 111 : finally, it is not confined in the x direction. The 
system of width A along the x direction is shown in Fig. 1, 
with periodic boundary conditions imposed along the x-axis. 
Within the analytical models the ribbon may be treated as ini- 
tially uniform in the y direction with negligible influence of 
the trap edges. Still, finite length of the ribbon is required in 
the numerical solution to the Gross-Pitaevskii equations. 

In the case of zero external forcing, linear stability analysis 
lf36h predicts a perturbation mode, which may be interpreted 
as quantum capillary waves with the frequency £2 £ depending 
on the wave number k — 2%jX as 
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where tilde indicates dimensional values and «o is the con- 
densate density at the unperturbed infinitely thin interface (i.e. 
in the Thomas-Fermi approximation). Taking into account a 
time-dependent magnetic force pushing BECs to each other, 
we obtain an equation for the linear interface perturbations 
lIU 
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where £(f ) is amplitude of deviation of interface from hydro- 
static equilibrium, pis is the Bohr magneton and B'{t) indi- 
cates the external (Stern-Gerlach) force due to the magnetic 
field gradient. 

We will consider a harmonically oscillating force, which 
may lead to the parametric instability. For simplicity of des- 
ignations, we use the standard scalings for a trapped system 
with the length scales measured in units of the oscillator char- 
acteristic length fl y = ^h/mcOy, time in units of (2(O y )~ l , and 
wave functions in units of y^o". Within these designations, the 
dimensionless system size is Ro = ^J2gh()/Ti(0y and the dis- 
persion relation for the capillary waves is £2 2 . = ^/yRolc' . In 
the same manner we define the dimensionless external force 
b(t) = iiBB'(t)a y /TK£>y, which we take in the form 




FIG. 1: (Color online) The system geometry: density distribution for 
two immiscible BECs, n\ and nj, for Rq = 30, with initial pertur- 
bations at the interface. The BECs are tightly trapped in z direction 
(means quasi-2D geometry), trapped with resulting TF profile in y 
direction, and is uniform in x direction. Dimensionless units are de- 
fined in Sec II. 



The approximate analytical solution to Eq. (O with the oscil- 
lating force may be obtained using the method of [37], which 
describes exponential growth of interface perturbations as 



£ = £osin(i2f/2 + (p)exp(af) 



(4) 



with the growth rate a, some amplitude £o, and the phase shift 
<p with respect to the driving force. Mark that the growth of the 
parametric instability is accompanied by interface oscillations 
with frequency £2/2. Then the perturbation growth rate may 
be found analytically using the equation 
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see [20] for details; the growth rate is plotted vs the perturba- 
tion wave number in Fig. 2 for £2 = 8.55 and different am- 
plitudes of the driving force. As we can see, maximal growth 
rate is achieved for the frequency of the driving force equal 
doubled capillary frequency, £2 = 2£2 £ , which is the condition 
of a parametric resonance. In that case the growth rate is given 
by a simple formula 
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b{t) = & c sin(£2*). 



(3) 



and the parametric instability may be excited even by an ex- 
tremely weak force. For other frequencies of the driving force, 
out of the resonance £2 ^ 2£2 £ , a finite force amplitude b c is 
required to produce the instability. We point out that in the 
case of moderate amplitude of the driving force, the instabil- 
ity domain in wave numbers is strongly localized around the 
resonance point, see Fig. 2. Another interesting character- 
istic of the linear stage is the relative value of the instability 
growth rate a and the perturbation frequency £2/2 in Eq. (0|, 
which shows how fast the perturbations grow in one oscilla- 
tion period of the driving force. As we will see below, the 
instability parameter 2a/£2 determines strength of the pro- 
cess and qualitative regime of the parametric resonance at the 
nonlinear stage. 

Let us qualitatively discuss dynamics of the interface ex- 
pected at the nonlinear stage of the parametric resonance, for 
the cases of weak instability (2a/£2 -C 1) and strong instabil- 
ity (2a/£2 3> 1). The transition between these two regimes 
happens for 2a/£2 ~ 1. In the case of weak instability one 
should expect nonlinear stabilization of the parametric reso- 
nance because of nonlinear modifications of the oscillation 



FIG. 2: Scaled growth rate of the parametric instability vs the wave 
number of the interface perturbation between BECs for £2 = 8.55 and 
different amplitudes of the driving force, obtained in 1I20I1 . Dimen- 
sionless units are defined in Sec II. 



FIG. 3: (Color online) Density distribution of BEC 1, n\, for/?o = 30, 
b c = 3 at t = 2.0, showing the maximal interface deformation in the 
regime of non-linearly stabilized parametric resonance. Dimension- 
less units are defined in Sec II. 



frequency B37I1 . Taking into account small but finite ampli- 
tude of the oscillations £0, one finds the oscillation frequency 
modified from Q./2 to Q./2 — kQ, where K is some factor, 
which may be found, e.g., from the numerical solution. Since 
the frequency changes together with the growth of the pertur- 
bation amplitude, then the system eventually leaves the insta- 
bility domain and the perturbation growth is saturated. Sub- 
stituting the modified frequency into Eq. (0 and taking a = 
we evaluate the maximal oscillation amplitude attained in the 
resonance as 



Co = y/b c k/2£lK. 



(7) 




In the opposite limit of strong instability, 2a/D.^$> 1, the ap- 
proach of small corrections to the linear solution employed in 
113711 fails. In that case we expect considerable growth of the 
perturbation amplitude already on one oscillation period with 
strongly nonlinear effects resembling qualitatively the RT in- 
stability, though with time-dependent effective gravity. 



in. NUMERICAL SOLUTION 

We investigate the parametric resonance numerically by 
solving the coupled Gross-Pitaevskii equations 
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where A r-V = d 2 x + d 2 Y , and the potentials Vi^{t,x,y) in- 
clude both the trapping potential and the driving potential 
as Vj = mco 2 (x 2 +y 2 ) /2 + (-iy'^ B B'(t)y/2 = V, + {-iy'V d . 
The trapping potential along z-axis is supposed to be much 
stronger than the two-particle interaction energy, and this 
leads to renormalization of the interaction parameters rf35ll . In 
all calculations, we keep the scaled system size Rq = 30 and 



FIG. 4: (Color online) Amplitude £ (t) of the interface deviation from 
equilibrium for Rq = 30 and b c = 1, 3 and the scaled amplitude of the 
driving force (dots). Dimensionless units are defined in Sec II. 



the repulsion parameter 7 = 0.01, which implies a sufficiently 
small thickness of the interface ~ 1/ y/fRo ~ 0.3 and a much 
smaller healing length - 1/R Q rj 3.3 • 10~ 2 Hfl]. 

The unper- 
turbed planar distribution of BEC density has been obtained 
numerically similar to |[TT[[T9ll . 

At t = we impose a cosine-shaped perturbation at the in- 
terface as shown in Fig. 1, which corresponds to the capillary 
mode with the perturbation wave length A = 3.44, and im- 
mediately turn on the driving magnetic (Stem-Gerlach) force, 
and study dynamics of the system for different amplitudes 
of the force b c = 1 — 10. Since we are interested in the 
parametric resonance, then we chose the external force fre- 
quency i2 = 8.55 equal to 2Q. C as calculated for the wave- 
length A = 3.44 of the initial perturbations. In order to define 
position of the interface in the study, we chose the density 
level equal to the density of the unperturbed solution at z = 0. 

System dynamics at the nonlinear stage of the resonance 
depends qualitatively on the magnitude of the external force 
b c - In the case of relatively weak force, e.g. for b c . = 1, 3, 
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FIG. 5: (Color online) Phase diagram of vortex generation in the 
fixed-size (fixed Rq) system, showing possible regimes with well- 
defined interface. For b c > 10 the system is not well-described in the 
mean-field approximation because it acquires too much incoherent 
excitations. Mean deviation of the interface from equilibrium £n, av- 
eraged on twice as many periods of modulations, as shown in Fig. 4, 
for various b c (time-averaged for irregular dependence) Rq = 30. The 
solid line corresponds to the analytical model (0 with K = 0.96. The 
dashed line separates two qualitatively different nonlinear regimes. 
Dimensionless units are defined in Sec II. 



the shape of the interface resembles the initial cosine-shaped 
perturbations, though with increased amplitude. As an exam- 
ple, maximal distortion of the interface for b c = 3 is demon- 
strated in Fig. 3 for the time instant t = 2.0. Time dependence 
of the perturbation amplitude for b c — 1,3 may be character- 
ized as rather regular modulated oscillations, see Fig. 4, with 
frequency ps £2/2 and with the maximal interface distortion 
attained after 3-4 oscillations. After that, dynamics of the in- 
terface deviation amplitude is reproduced periodically in the 
form of modulations. We point out that the force magnitudes 
b c = 1, 3 belong to the limit of weak instability with the insta- 
bility parameter 2oc/£2 = 0.095, 0.29, respectively. Thus, in 
agreement with the expectations of the analytical model, the 
parametric resonance is stabilized at the nonlinear stage in the 
case of weak instability. Numerical simulations indicate that 
the oscillation amplitude increases with increasing magnitude 
of the driving force as shown in Fig. 5: comparison of the 
numerical results with Eq. (|7) demonstrates good agreement 
and allows evaluating the factor K for the chosen problem pa- 
rameters as K — 0.96. 

Taking larger force magnitude b c = 6 we obtain the instabil- 
ity parameter comparable to unity, 2a/£2 = 0.57, which leads 
to considerable nonlinear effects in the interface dynamics 
presented in Fig. 6 for selected time instants t = 0.08 — 3.46. 
Starting with small initial perturbation, the interface oscillates 
with noticeably smaller frequency than £2/2. Besides, the os- 
cillations are accompanied by strong broadening of the inter- 
face on the reversal motion of humps transformed into hol- 
lows, e.g. at t = 1.14. The second growth of the hump on 
the snapshots t = 1.66 — 2.08 is characterized by a compli- 
cated shape of the interface, quite distinct from the original 
harmonic perturbation. We mark resemblance of the interface 
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FIG. 6: (Color online) (Above) The scaled amplitude of the driving 
force, and amplitude of the interface deviation from equilibrium for 
Rq = 30 and b c = 6; the markers indicate the time instants for the 
density snapshots below. (Below) Snapshots of density distribution 
of BEC 1 n\. Dimensionless units are defined in Sec II. 



shape at t = 2.08 and the "mushroom" structures for the RT 
instability 1 181 |20ll . Due to the large amplitudes we observe 
the tendency of the mushroom cap to detach from the main 
bulk of the BEC at t = 2.4 and to form a large blob, i.e. a 2D 
soliton, according to the mechanism discussed in 111 911 . The 
mechanism of the blob detachment is also related to the cap- 
illary instability studied in ||2lll . Collapse of the mushroom 
structure may be also compared to the parametric resonance 
observed in the combustion system of a flame front in the ef- 
fective gravitational field produced by sound ll33l l34ll . Still, 
in the case of b c = 6, detachment of the blob is not complete, 
and it is followed by reattachment at t = 2.66. Because of the 
detachment-reattachment process, the definition of the oscil- 
lation amplitude becomes ambiguous, which leads, in partic- 
ular, to irregularities in the plot of £(f). Reattachment of the 
blob to the bulk of BEC is accompanied by momentum trans- 
fer, i.e. bounce, which produces a grey soliton traveling in the 



FIG. 7: (Color online) Amplitude of the interface deviation from 
equilibrium for Rq = 30 and b c = 8 — 10, and the scaled amplitude 
of the driving force. The time instants of vortex detachment are indi- 
cated by circles. Dimensionless units are defined in Sec II. 



bulk of BEC to the left in the snapshot t = 3.46. The soliton is 
formed because mean squared deviation of the external force 
is well above the critical value b cr = yRo /2 discussed in detail 
in ifTTIl . At the same time, interface dynamics for b c — 6 does 
not lead to production of vortices and most of the effects, ex- 
cept production of grey solitons, may be described within the 
quasi-classical approach. 

Generation of vortices may be observed at larger magni- 
tude of the driving force, b L ■ = 8 — 10 in our simulations. The 
approximate boundary between two regimes is indicated in 
Fig. 5 by the dashed vertical line, and it corresponds to the 
instability parameter about 2a /SI = 0.8, i.e. close to unity, in 
agreement with the expectation of the analytical model. Fig- 
ure 7 presents amplitude of the interface deviation from equi- 
librium for b t = 8 — 10 until the instants of droplet detachment 
and formation of quantum vortices in BEC 1 with the vortex 
core filled by BEC 2 and vice versa - skyrmions. After that 
instant, definition of the interface deviation becomes ambigu- 
ous. As we can see in Fig. 7, skyrmions are generated after 
several interface oscillations for b c . = 8, 9, and on the very 
first oscillation for b c = 10. Similar to Fig. 6, the oscillations 
are extremely irregular with characteristic frequency notice- 
ably smaller than SI/ 2. Taking the simulation run for b c = 8 
as an example, we observe the front shape in the oscillations 
for t < 4.5 resembling qualitatively the structures presented in 
Fig. 6. For this reason, we focus on the snapshots for b c = 8 
and t > 4.5 in Fig. 8 when qualitatively new quantum effects 
come to play. In particular, Fig. 8 demonstrates density and 
phase distributions of BEC 1 for the time instants t = 4.8 illus- 
trating detachment of the first droplets (2D solitons); t = 4.94 
when the soliton splits into the vortex-antivortexpair; t = 5. 12 
when the vortex- antivortex get almost coalescent because of 
mutual attraction; t = 7.7 when the vortex-antivortex are sep- 
arated again and drifted away from the main bulk of BEC2. 
In addition, new vortex-antivortex pairs may be observed at 
f = 7.7. 
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FIG. 8: (Color online) Snapshots of density distribution of BEC 1, 
n\, and phase distribution for Rq = 30 and b c = 8, showing vortex 
generation from the interface between BECs. The detached droplet 
on (b) splits into a vortex-antivortex pair on (c), and this happens pe- 
riodically on (d,e) and further, due to periodic external force. Phase 
plots show topological excitations. Dimensionless units are defined 
in Sec II. 



IV. GENERATION OF VORTICES IN THE PARAMETRIC 
INSTABILITY 

A. Pseudo-spin description of the two-component BEC 

The regime of large-amplitude force is characterized by 
generation of quantum vortex-antivortex pairs, so that each of 
the vortices carries circulation, while total circulation is zero. 
In the 3D case, similar vortices would be ring-shaped, though 
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FIG. 9: (Color online) Cut of n\ and showing dynamic AT- vortex 
and antivortex of the 1 component (a 2D skyrmion pair) in Fig. 8(d) 
at t=7.7 and z=-9.635. The 2 component fills the vortex cores of the 
1 component. Dimensionless units are defined in Sec II. 



the Kelvin waves can also split the vortex lines into smaller 
pieces. The ground state and the vortex structure in the two- 
component BEC is determined by the interaction parameters. 
Since the system in study is characterized by weak repulsion, 
7= gn/g — 1 <C 1, the topological excitations are vortices in 
one component with the core filled by atoms of another com- 
ponent; such a structure is known as an Anderson-Toulouse 
(AT) vortex. The dynamic states shown in Fig. 8 (c, e) present 
clear examples of the AT vortex pairs; Fig. 9 shows a cut 
of BEC densities m 2 taken from Fig. 8 (e) at y = —9.635, 
along the line crossing the vortex axes. The obtained vortices 
may be also described as skyrmions, though the term skyrmion 
refers to a rather wide class of topological solitons. Histori- 
cally this term originates from the theory by T. H. R. Skyrme 
[38], which is quite similar mathematically to description of 
multicomponent BECs ll39ll . For a small repulsion parameter 
7< 1, radius of the AT vortex core is given by §q / \/7> wnere 
<jjo = h/ \J2mgn, i.e. it is 1 / yfy times larger than the core ra- 
dius in a single-component BEC. This value is consistent with 
the vortex structure in Fig. 9, and it is also comparable to the 
effective width of the interface between the components rf20Tl . 

Let us describe generation of vorticity in our system. A nat- 
ural way to reveal dynamical difference between single- and 
multi-component BECs is to go over from the representation 
of the GP Lagrangian and the wave functions (yf\, y/2) to 
the representation of a normalized spinor 

s=ki«p(i/J/2), jfe«p(-i/J/2)] (10) 

and the variables ^/prexp(ia) corresponding to total density 
Pt and total phase a. The normalization implies x\ + x\ = 1> 
and therefore we can choose %\ = cos (% /2), %2 — sm Gk/2); 
the new variables p T ,a,j5,X ^ real-valued functions of 
coordinates and time. If the wave functions are presented 
as y/j = \/"7 ex P(^')' tnen P T = ni + ra 2. a = (0i + <fa)/2, 
J3 = 01 - <p 2 , and % = 2 arctan y/n 2 /n l . 

It is also convenient to define the pseudo-spin vector as S = 
s<7s, where a — [a x a y o z ) T , and o xy , are the Pauli matri- 
ces. Then one readily obtains S x = cos /3 sin^, S y = sin /3 sin^, 



S z = cos x, with the relation S x + S? + S 2 = 1 . The superfluid 
mass current pr^eff is derived from equations of motion, Eqs. 
©,©, as 



pTVeff = |Vi| 2 Varg(v/i) + |i//2| 2 Varg(v/ 2 ) 
= Pt 



va + —cosx 



- Pt 



Va + 5. 



s y vs x - s x vs y 
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see M40I1 for derivation of the last term. Equation ([Til demon- 
strates that vorticity V x \ e ft of the effective superfluid veloc- 
ity can be non-zero in a multicomponent system without sin- 
gular regions of the order parameter, in contrast to the single- 
component case. Using Eq. (fTTT i we compute evolution of 
vorticity in our system, and plot V x v e // in Fig. 10 for the 
same moments of time and the same parameters as shown in 
Fig. 8. This function may be also interpreted as density of the 
topological charge, and when integrated over space it gives 
the topological invariant. This invariant, or Pontryagin index, 
classifies the stationary states of multi-component BECs; in 
our system it is formed from the z-component of the vorticity 
pseudo-vector, which can be derived from Eq. (fTH as ll4Tll 

q(r) = ^e ,; S [diS x djS] = ^(Vx v eff ) z . (12) 

In a more general case the proof of Eq. (fT2t can be obtained 
by considering the form XnXh, where a, b count the spinor de- 
grees of freedom 04111 . Topological charge of a single local- 
ized AT vortex is Q — J d 2 rq(r) = 1. 

One observes that vorticity V x v e // is generated on the in- 
terface between BECs. In our superfluid system the interface 
region may behave locally like a rotating solid body during 
some finite time intervals. In particular, such behavior may 
be observed in Fig. 10 (a) for t = 4.4 , which precedes the 
detachment of two vortex- antivortex pairs, Fig. 10 (b). When 
the pairs detach, the interface acquires an opposite effective 
vorticity sign, Fig. 10 (c), and then, after a period of cap- 
illary oscillations, another pair of vortex-antivortex is gener- 
ated, Fig. 10 (d). In Fig. 10 we also observe excitation of 
compressible modes, which may lead to dissipations due to 
two-body processes. However, the dissipation processes are 
beyond the scope of the GP model and are not considered in 
the present work. 



B. Dynamics of vortices 

The vortex-antivortex pairs are formed from the bubbles, 
and propagate towards the BEC edge. The pairs are generated 
on the interface one after another, and they annihilate close 
to the BEC edge after some sporadic dynamics. In this sub- 
section we investigate dynamics of the first pair on the stage 
of propagation towards the edge. In our simulations interac- 
tion between different vortex pairs is negligible, and a typical 
trajectory of a pair is represented in Fig. 11, where we glue 
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FIG. 10: (Color online) Snapshots of topological charge density of 
the two-component BEC for Rq = 30 and b c = 8, for time instants 
the same as in Fig. 8. Topological charge of the system is generated 
on the interface between BECs. The direction of the z-axis is along 
x x y. Dimensionless units are defined in Sec II. 



together cuts from the full density profiles taken at 10 differ- 
ent time moments. The snapshots are limited by the dashed 
borders which show local density. 

First let us explain how vortex pairs are formed from the 
bubbles. Because of the non-linearity of the process of bubble 
detachment from the interface it is not clear whether one could 
predict the initial velocity of the bubble vo immediately after 
the detachment. However, instead we can predict the critical 
velocity of the bubble, above which it splits into vortex pairs. 
This velocity results from topology of the order parameter, be- 
cause quantized vortices are topological excitations. In the TF 
approximation, the bubble of the BEC 2 produces a cavity in 
the bulk of the BEC 1, which travels with the bubble velocity, 
and given that the bulk of BEC 1 is at rest, then, up to loga- 
rithmic accuracy, its phase must have a "jump" (discontinuity) 
at the cavity center. This jump can exceed 2n, and therefore 
give rise to formation of the vortices, only when diameter L 
of the bubble exceeds 2nh/mvQ. Therefore the critical bubble 
velocity of the BEC 2 above which the bubble splits into a 
vortex-antivortex pair, is 



2%h 



(13) 



where m\ is mass of atom of the BEC 1 (in our numerical sim- 
ulations we take m\ = = m). Equation (fT3l is valid both 
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FIG. 1 1 : (Color online) Snapshots of location of the first generated 
vortex-antivortex pair in the 1 component (at y < 0), with correspond- 
ing fillings of the cores of the first vortex-antivortex pair of the 2 
component (at y > 0), for Ro = 30 and b c = 8. Snapshots are taken 
at time instants (a) t = 5.0, (b) 5.4, (c) 5.8, (d) 6.2, (e) 6.6, (f) 7.0, 
(g) 7.4, (h) 7.8, (i) 8.2, and (j) 8.6, and only the area close the pair is 
shown for each pair (marked by dashed borders). Note that the time 
instant t = 5.12 from Fig. 8(d) is not shown in order to make the 
picture easier for understanding. The second pairs, generated after 
the first ones, are not shown for the same reason. At the time in- 
stant t = 8.6, shown on (j), dynamics of vortices is sporadic due to 
localization of the system and strong influence of the external field. 
Dimensionless units are defined in Sec II. 



for 2D and 3D geometries. We point out that our argument, 
leading to Eq. ( fT3l , is based on topology of the order param- 
eter and therefore is general. The only assumption made is 
validity of the TF approximation L > 2Ai„ t . The result Eq. 
(fob can be applied as well for a single component BEC when 
it is stirred by an obstacle with the size much larger than the 
healing length of the BEC. As a result, the number of pairs 
generated by a bubble with a given initial velocity v>o is equal 
to the integer part of vo/v cr . In our system for^o = 30, b c = 8 
we observe formation of a single pair from each bubble. We 
measure velocity of the first pair as w 1 .25 ^licOy./m, which 
is slightly higher than predictions of Eq. (fT3~t computed with 
L = 2A,„, as v cr ps Q.95 yJhcOy/m. In general one should take 
into account compressibility of the bubble B42I1 ; the capillary 
instability [21] for large bubbles may influence the results too. 

Equation ([T3l can be used to demonstrate that the lowest 
critical velocity that breaks the Landau criterion of superflu- 
idity, is related to topological excitations rather than to sound- 
like excitations. Indeed, the sound speed c s in the trap cen- 
ter is equal to V2Rq, in the dimensionless variables speci- 
fied in Sec. II. It should be compared with the dimensionless 
Eq. ( fT3l l, 4n/L, where L is size of the obstacle in units of 
a,.. Using L = 2/(-y/yR ) with R = 30, and y = 0.01, rele- 
vant for droplets seen in evolution of our system, we obtain 
v cr 0.44c 5 , and this is confirmed by the observation that the 
detached droplets split into vortex-antivortex pairs. 

Now we turn to the trajectories of the pairs, which are es- 
sentially straight lines modulated by relatively weak oscilla- 
tions. The straight lines indicate the stationary propagation 
of the pair. Indeed, an isolated stationary vortex pair (vor- 
tex ring in 3D) moves with a constant velocity with respect to 
the quiescent bulk when the attraction Magnus force between 
the vortices in the pair is compensated by the repulsion Mag- 
nus force due to the vortex motion. The net Magnus force is 
equal zero when velocity of the vortex pair is the same as the 
velocity created by one vortex of the pair in the position of 
the other. The system is invariant with respect to the Galilean 
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coordinate transformation as can be readily shown in the GP 
model. However when the atomic cloud, containing a vortex, 
is subjected to acceleration, the vortex exhibits inertia due to 
mass of the core filling 



with e 
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-e Jl , i,j = x,y, and Vi = ±1 is the wind- 



nif = mjn%ts. i} 



(14) 



where A inf = Ro/^/Y is the radius of the vortex core, and we 
have chosen to consider a vortex in the 1 component filled by 
the 2 component. On the other hand, the effective mass of the 
vortex itself may be obtained taking into account the relativis- 
tic character of the collective excitations of BEC coupled to 
the vortex state HH: 



(15) 



where E v = nnih 2 \n (R tot /Ai nt ) jm\ is energy of a vortex in 
the BEC 1 in a cylindrical sys tem lim ited by a container with 
a large radius R tot , and c s — y^gn/m is the sound speed in the 
trap center. Thus the effective action, describing the degree 
of freedom of the system associated with the position of the 
vortex line element R(f ), contains the "kinetic" contribution 
of the form 



Ski, 



dt (m v + m/)R 2 /2. 



(16) 



The "potential" contribution describes forces acting on a vor- 
tex. It includes the Stern-Gerlach force applied to the filling 
composed of atoms of the 2 component (in our system the 
force is directed along y) 



dt^Ry^. 

mi 2 



(17) 



and the lift force, i.e. the quantum analogue of the Mag- 
nus force, which may be derived from the GP model. In the 
derivation we start from the action for a two-component BEC 



S= dt 



d 3 r 



.J=1.2 



(18) 



The "potential" contribution to the action associated with in- 
teraction of the vortex with the background velocity is 

S M = f Ac/ 2 rRe[i>r(r-R(f))<3 f V/i(r-R(f))] = 



(19) 



dtd 2 r [^*(r-R(f))V^i(r-R(f))- 



V / 1 (r-R(r))V V / 1 *(r-R(f))]R, 



where y/i is the stationary wave function of a vortex (the 
chemical potential was dropped for brevity because it does not 
alter the result). The wave function Yi(r — R(f)) is Taylor- 
expanded for small R(f ) as 

V/i(r-R(0) = V /i(r)+R(f)V V / 1 (r) + 0(R 2 ,...), (20) 

and only the first order in R(f ) is kept. After that the spatial in- 
tegration is readily done, and we obtain for a singly-quantized 
vortex in the 1 component 



Sm — —V\Tthn\ J dte' j RiRj, 



(21) 



ing number of the vortex with the effective vorticity along z 
(+1), or in the opposite direction (—1). It follows from the 
energetic conditions, and from our real-time numerical simu- 
lations, that vortices with |vi | > 1 are not formed during the 
evolution of our system. Thus, vortices and anti-vortices ex- 
perience opposite forces due to a background superfluid flow 
with velocity R. As a result, a single AT vortex in a homoge- 
nous two-component BEC is described in 2D by the action 

5 V [R,R] = S ext + j dt [{m v + m f )R 2 /2 - Vinhne'-'RiRj] . 

(22) 

We note that the action of the external force on the 1 BEC is 
not explicit in Eq. d22b . because it implicitly contributes to R 
by inducing collective bulk oscillations. 

Equation (l22l describes the quantum analogue of the Mag- 
nus force, which acts perpendicularly to velocity of a vortex 
R(f). Using Eq. d22"i i and going to dimensionless variables 
specified in Sec. II, we obtain the dimensionless form of the 
force acting on a vortex in the 1 BEC with non-singular core 
filled by the 2 BEC 



d 2 R, 



4Vijgo£% + fcKO 
' l+2yln(yJW4o) ' 



(23) 



where bi(t) are Cartesian components of the external force 
applied to the 2 BEC, b y (t) = b(t), see Eq. ©. Equation ( 1231 



shows that for systems with size R, ot < (2/ y jRq) exp(l /2y) 

the vortex mass is dominated by mass of the core filling. In 
our system R tot ps Rq, and the vortex mass is approximately 
that of the filling. We conclude that the oscillations of the pair 
size observed in Fig. 11, are due to the bulk motion driven 
by the external force, and due to the external force applied to 
the filling. Motion of the bulk determines bulk velocity with 
respect to vortices, and hence the Magnus force. As follows 
from Eq. d23l . the frequency of the oscillations is close to 
the frequency of the external force, i.e. ps Q. = 8.55, which is 
confirmed by the numerical solution shown in Fig. 1 1 . 



V. SUMMARY 

We have studied the parametric instability at the interface of 
two immiscible BEC components pushed towards each other 
by an oscillating force. The instability develops due to the 
parametric resonance, which pumps quantum capillary waves 
at the interface. At moderate amplitudes of the driving force 
the instability is stabilized at the nonlinear stage due to mod- 
ifications of the doubled oscillation frequency in comparison 
with frequency of the driving force. In that case the BEC in- 
terface demonstrates oscillations with modulated amplitude, 
which depends on the strength of the driving force. When the 
amplitude of the driving force is large enough, we observe de- 
tachment of droplets from BEC components and generation of 
quantum vortices - skyrmions. The skyrmions are born as vor- 
tex pairs and move almost stationary from the interface to the 
trap edge. The properties and dynamics of the skyrmion pairs 
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are discussed. We have introduced the critical velocity v cr of a 
droplet of one BEC in a bulk of another BEC from topological 
arguments, and derived the analytical formula for the quantum 
counterpart of the Magnus force acting on skyrmions. 

Note added in the proof. Recently, the authors became 
aware of the work [45] on the Faraday resonance in two- 
component BEC confined to a quasi ID geometry. The main 
difference between the current work and Ref. [45 ] is in the ge- 
ometry (multidimensional versus ID, respectively) and in the 
nature of the surface waves amplified by the oscillating exter- 
nal force. Here we study interfacial capillary modes on the 
common surface of two BECs, while in the work Il45ll these 



are the surface modes at the edges of each BEC facing vac- 
uum. 
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